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Introduction 

The successful design of advanced com- 
posite aerospace structures requires a detail 
knowledge of strains and stresses that are in- 
duced by the applied loads. The inherent 
strength weakness between adjacent plies con- 
stitutes a weak link in a laminate which under 
the action of interlaminar stresses may result 
in the dominant interlaminar failure mode 
know as delamination. The delamination may 
then yield the ultimate laminate failure. Con- 
sequently, predicting accurate interlaminar 
strains and stresses is an essential facet of 
composites analysis. An comprehensive ac- 
count of the currently available methods for 
the interlaminar stress recovery is presented 
by Byun and Kapania. 1 

The first-order shear deformation theory 
with the use of constitutive relations produces 
highly inaccurate interlaminar stresses that 
exhibit nonphysical discontinuity at the ply 
interfaces of composite laminates. This diffi- 
culty is overcome by integrating three-dimen- 
sional equilibrium equations of elasticity the- 
ory to recover the interlaminar stresses. 2 This 
alternative procedure, which is not entirely 
consistent with the assumptions of the first-or- 
der shear deformation theory, provides rea- 
sonably accurate predictions for the interlami- 
nar stresses for closed form analytic solutions. 
When the procedure is applied to the results of 
finite element analysis, the accuracy of the 
interlaminar stress predictions tends to be 
rather poor. The main reason for this is that 
the integration of three-dimensional equilib- 
rium equations involves expressions of 


strain gradients whose recovery from the finite 
element shape functions is known to be infe- 
rior. 

This report discusses an accurate recovery 
procedure of interlaminar shear stresses from 
finite element solutions to laminated sandwich 
and composite plates. The procedure is gen- 
eral enough to be useful as a full-field post- 
processor of all finite element response quan- 
tities. 

Approach 

The preferred procedure for computing 
through-the-thickness distributions of inter- 
laminar (transverse) stresses in laminated 
plates and shells is by means of integrating 
three-dimensional equilibrium equations of 
elasticity theory. The integrals contain partial 
derivatives of strain measures of plate/shell 
theory. These quantities are independent of the 
thickness coordinate and are commonly at- 
tributed to a reference plate/shell surface. For 
example, upon integration of the elasticity 
equilibrium equations, the interlaminar shear 
stress, , corresponding to the first-order 
shear-deformation plate theory can be written 
as 

**=Qr«* + Q«-*., 0) 

where a comma stands for partial differentia- 
tion; e is the vector of strain measures, 

e = {£„, e yo , Yxyo . K xo . K yo - < xyo 1 ( 1 a ) 



Q j are the vectors whose components are 
functions of the thickness coordinate, z, 

Q j = (lb) 

where j = 1,6 and 

= -jcf Py = -jC-j k) zdz (lc) 

with cjj k) denoting the elastic stiffness coeffi- 
cients of the k-th ply. 

The difficulties in obtaining reliable pre- 
dictions of strain gradients from a finite ele- 
ment analysis are well established. For exam- 
ple, in first-order shear deformation elements 
based on Reissner-Mindlin type theories, the 
displacement (kinematic) fields are approxi- 
mated with C° -continuous shape functions en- 
suring continuity of displacements along ele- 
ment boundaries. The strains are represented 
by the first gradients of displacements; hence, 
they are only C ' 1 continuous, exhibiting non- 
physical discontinuities along element bound- 
aries. Only at a limited number of discrete lo- 
cations within the interior of a finite element 
reasonably accurate strain results can be calcu- 
lated, these are referred to as optimal points. 
The first gradients of strain are obviously even 
less accurate, and hence reliable strain deriva- 
tives cannot be obtained directly from the el- 
ement interpolation functions. 

Tessler and co-workers 3 ' 6 recently devel- 
oped a smoothing variational formulation 
which combines discrete least-squares and 
penalty-constraint functionals in a single 
variational form. The strain (or stress) compo- 
nent to be smoothed and its Cartesian gradi- 
ents enter in the variational principle as inde- 
pendent variables. For a sufficiently large 
value of the penalty number, the approach en- 
ables the resulting smooth strain field to be 
practically C 1 -continuous throughout the do- 
main of smoothing, exhibiting superconver- 
gent properties of the smoothed quantity. The 
continuous strain gradients are also obtained 


directly from the solution. The use of the re- 
covered stress (strain) gradients has been 
demonstrated for assessing errors in finite el- 
ement analysis by computing the residual in 
the strong form of the equilibrium equation. 
Also, the method has the versatility of being 
applicable to the analysis of rather general and 
complex structures built of distinct compo- 
nents and materials, such as found in aircraft 
design. For these types of structures, the 
smoothing is achieved with “patches”, each 
patch covering the domain in which the 
smoothed quantity is physically continuous. 
The mathematical basis of the smoothing 
methodology is briefly summarized below. 

Let e p =e(x p ) ( x p =(x' p ,xp , p=l,2,...,N) 
represent a set of discrete optimal values for a 
specific strain component defined on the two- 
dimensional domain f2 = [x = (^',^ 2 )c5t 2 ], 
with these strains resulting from a finite ele- 
ment analysis. We seek a smooth, continuous 
strain e(x) and its gradients e, =-^e (i=l,2) 
which can best represent the discrete data. 

The problem can be cast as a minimization 
of a least-squares/penalty-constraint func- 
tional, such that 

80 = 0 

where 0 is given by 

**•^1 [Cp-eCx ,,)] 2 

P= ' ( 2 ) 
+Af [(e l ~e i ) 2 +(e 2 -e 2 ) 2 ]da 

a 

where N is the total number of data points, 
x p denotes the position vector of the data 
point, A is a penalty number, and 0, and 0 2 
are independent continuous functions also de- 
fined on £2. The minimization of 0 is per- 
formed with respect to the coefficients in e 
and 0, which serve as the unknowns. The first 
term in Eq. (2) represents a discrete least- 
squares functional with the term 
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A p =£ p -£(*p) 


(3) 

signifying the error in e p as compared with 
the smooth solution e(x p ). The second term in 
Eq. (2) is a penalty constraint functional 
which, for sufficiently large values of A , 
yields the following constraint relations 

e, -» 6, (/ = 1,2) on £2 (4) 

The proper interpretation of Eq. (4) is that for 
all practical purposes the 6, (i = 1,2) represent 
the gradients of e. The variational principle 
in Eq. (2) is readily discretized with (^-con- 
tinuous element approximations, forming a 
mesh of “smoothing” elements. The smooth- 
ing finite element analysis then provides con- 
tinuous and reliable solutions for the strain 
and its gradients. The latter quantities are then 
employed in Eq. (1) for the recovery of inter- 
laminar shear stress distributions through the 
thickness. 

In this paper, the accuracy and robustness 
of the two-dimensional smoothing discretiza- 
tion is examined for the problem of recovering 
accurate strain gradients, and subsequent inte- 
gration of equilibrium equations to compute 
interlaminar shear stress distributions in lami- 
nated composites. The model problem is a 
simply-supported rectangular plate under a 
doubly sinusoidal transverse pressure load. 
The problem has an exact analytic solution 
which serves as a measure of goodness of the 
recovered interlaminar shear stresses. 

The following charts describe the essential 
theoretical features of the approach and sum- 
marize the numerical results. 
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Eight Strain Measures: 

First Partial Derivatives of Kinematic Variables 




e y,x 

Ko 

► = * 

Q x,y 

y 

—o J 


u x,x ~ u y,y 


W .x + 9y 


w<y + e x 


Inplane Strains 


Bending Curvatures 


Transverse Shear Strains 
(uniform through thickness) 


FEM: C 0_ continuous interpolations u , v f w, 0 xt 0 y 


Integration Scheme 


(i) 

= Ql e ,x + Q6 ' e ,v 

Shear Stress 

(la) 

® { ^xo » ^yo * Y xyo 1 ^xo * * ^*xvo ) 

Strain Vector 

(1b) 

Q; “ (*7l ; * Qlj * ^6 J * Pi j* Pi j ' P6j} 

Material-Thickness 
Coordinate Vector 

(1c) 

% =-\A k)dz ’ Pij=-l z ^ )zdz 

Qj components 
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Smoothing Analysis 


(i) i.=i(x.) 


(2)0=±Z[S p -e(x p )f 

™ p= 1 


Sample optimal FEM 
discrete strain values at Gauss 
points of finite elements 


+AJ l(e A -e i ) 2 +(e 2 -e 2 ) 2 ]di2 


(3) (/ = 1,2) 


Apply Least-squares/ 
penalty-constraint finite 
element method to smooth 
FEM strains for entire domain 


Recover C°-continuous, 
smoothed strain gradients 
throughout FEM domain 


Numerical Examples 

SS Square Plate Under Doubly Symmetric 
Sinusoidal Normal Pressure 


4 

y 

i 

SS 


• Loading: q sin(nx/L) sin(ny/L) 




• Span/Thickness=30, 10, and 5 

SS 

L q(*,y) 

L 

SS 

Lamination 


SS 

^ X 

• Gr/Ep: [0/90/0] s 6-ply 

• Gr/Ep: [O^OyOJs 18-ply 

• Sandwich: [(Q/90)s/hJ s 21-ply 
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FSDT vs. 3D Elasticity 
Interlaminar Shear Stresses 


Thin, 6 Ply Gr/Ep Laminate 
[0/90/0]$, L/2/i = 30 



FSDT : First-order Shear Deformation Theory Solution 

3D ELASTICITY: 3D Elasticity Theory Solution 


FSDT vs. 3D Elasticity 
Interlaminar Shear Stresses 

Moderately Thick, 1 8 Ply Gr/Ep Laminate 


[0 3 /90 3 /0 3 ]$, L/2/i = 10 
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FSDT vs. 3D Elasticity 
Interlaminar Shear Stresses 


Thick, 21 Ply Sandwich Laminate 
[(0/90) 5 /A c ] s , L/2h = 5 



FEM and Smoothing Meshes 


FINITE ELEMENT MESH 
9-NODE ANS SHELL ELEMENTS 



+ Optimal Curvatures @ 
2x2 Gauss Points 


SMOOTHING ELEMENT MESH 
3-NODE SMOOTHING ELEMENTS 



• K, 0j,0 2 DOF 
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Finite Element Results 

Thin, 6 Ply Gr/Ep Laminate 
[0/90/0]$, L/2/j = 30 


Response 

Quantitv 

% Error in 
Max. Value 

Deflection 

0.02 

Rotation 

0.03 

Curvature 

4.69 

Interlaminar Shear 

8.80 

Interlaminar Shear T yz 

8.20 


FEM vs. Recovered 
Interlaminar Shear Stresses 

Thin, 6 Ply Gr/Ep Laminate 
[0/90/0]$, L/2h = 30 
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FEM vs. Recovered 
Interlaminar Shear Stresses 

Thick, 21 Ply Sandwich Laminate 
[(0/90) 5 //j c ] 5 , L/2h = 5 



Results Significance 


• Direct evaluation of interlaminar shear 
stresses from FEM underestimates their value 
considerably (unconservative for design) 

• FEM-based smoothing method recovers 
improved interlaminar shear stresses that are 
slightly overestimated (conservative for 
design) 

• Recovered interlaminar shear stresses are 
improved by at least a factor of 3 [Error of 
2.8%(Recovery) vs. 8.2%(FEM)] 
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Summary 


• Strain smoothing produces high-accuracy 
strains and their gradients that are C°- 
continuous over full domain 

• ‘Smoothed’ strain gradients yield superior 
recovery of interlaminar shear stresses over 
direct FEM recovery 

• Effective for laminated composite and sandwich 
plates 

• Viable full-field post-processor of all strains & 
stresses in general purpose FE software 
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